Abstract 

In this paper we develop an analytical framework for the study of 
electrochemical impedance of mixed ionic and electronic conductors 
(MIEC). The framework is based on first-principles and it features 
the coupling of electrochemical reactions, surface transport and bulk 
transport processes. We utilize this work to analyze two-dimensional 
systems relevant for fuel cell science via finite element method (FEM). 
Alternate Current Impedance Spectroscopy (AC-IS or IS) of a ceria 
symmetric cell is simulated near equilibrium condition (zero bias) for a 
wide array of working conditions including variations of temperature 
and H2 partial pressure on a two-dimensional doped Ceria sample 
with patterned metal electrodes. The model shows agreement of IS 
curves with the experimental literature with the relative error on the 
impedance being consistently below 2%. Important two-dimensional 
effects such as the effects of thickness decrease and the influence of 
variable electronic and ionic diffusivities on the impedance spectra 
are also explored. 
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1 1 Introduction 

2 Mixed ionic and electronic conductors (MIEC), or in short mixed conductors, 

3 are substances capable of conducting both electrons and ions, and for that 

4 reason they are used in many applications, most notably in catalysis and 

5 eletrochemistry: they have been employed in gas sensors, fuel cells, oxygen 

6 permeation membranes, oxygen pumps and electrolyzers. 

7 The study of the alternate current properties of MIEC aides in under- 

8 standing many of the physical chemical phenomena related to the behavior 

9 of defects, electrochemistry and interfaces. A technique frequently used to 

10 probe the interplay between these processes is impedance spectroscopy (IS). 

11 IS consists in injecting a "small" sinusoidal current into an electrochemical 

12 sample, a fuel cell for example, which is initially under steady-state condi- 

13 tions. This perturbation in turn induces a small sinusoidal and de-phased 

14 perturbation of the voltage. From the measurements of voltage and current 

15 over a wide set of frequencies, one can compute the complex impedance of 

16 the system. When the experiment is compared against a suitable model, 

17 impedance spectroscopy helps understand the linear physics of electro-active 

18 materials. 

19 The tools used to deconvolute impedance spectra and relate them to 
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1 physical-chemical quantities are usually limited to one- dimensional equiva- 

2 lent circuits^.- Even though the ID approach is very useful because it enables 

3 the comparison of different processes, it sometimes fails to help interpret sat- 

4 isfactorily physical chemical phenomena that extend to several dimensions. 

5 Only a handful of works attempted to scale up to two dimensions, and gen- 

6 erally have been constrained to the steady-state setting— .- 

7 In this paper we develop a fast method for the computation of impedance 

8 spectra for highly-doped mixed conductors in a 2D setting under geometri- 

9 cally symmetric conditions. The system studied was chosen so that it is not 

10 too cumbersome algebraically and readily relatable to experiments. However 

11 the methodology is general and it can be easily extended to 3D, to dis- 

12 symmetric systems under non-zero bias and to complex chemical boundary 

13 conditions. 

14 The paper proceeds as follows: we first develop a model for impedance 

15 spectroscopy and determine the impedance equations,- then we compare our 

16 results to experimental data, finally we study the influence of parameter 

17 variation on the IS: the thickness of the sample, the rates of the chemical 

18 reactions at the exposed MIEC surface and the diffusivity profiles. 

19 After non-dimensionalization of the full drift diffusion equations, we find 
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1 that the ratio between the Debye length and the characteristic length scale of 

2 the material is remarkably large, hence we singularly perturb the governing 

3 equations and we deduce that electroneutrality is satisfied for a large portion 

4 of the sample. Then we apply a small sinusoidal perturbation to the potential, 

5 which mathematically translates into a regular perturbation of the equations; 

6 after formal algebraic manipulations we collect first order terms and deduce 

7 two complex and linear partial differential equations in 2D space and time. 

8 Thank to linearity, the Fourier transformation of these equations and their 

9 boundary conditions leads to the determination of the complex impedance 

10 spectroscopy equations which we solve in 2D space for the frequencies of 

11 interest . 

12 We verify our numerical results against experiments that are relevant for 

13 fuel cell applications. In particular, we study the case of a Samarium Doped 

14 Ceria (SDC) sample, immersed in a uniform atmosphere of argon, hydrogen 

15 and water vapor. The sample is symmetric and reversible and has been the 

16 subject of extensive research—.- We find excellent agreement between the 

17 computed impedance spectra and experimental data. This shows that the 

18 approximations and the model are likely to be valid, hence this framework 

19 could help address a number of important fundamental physical/chemical 
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1 issues in mixed conductors. 



2 2 System Under Study 

3 The physical system under study is a two-dimensional assembly which con- 

4 sists of a mixed oxygen ion and electron conductor slab of thickness 2/2 

5 sandwiched between two identical patterned metal current collectors, FigdJ 

6 The patterned collectors are repeated and symmetrical with respect to the 

7 center line Fi. Hence the system to be reduced to a repeating cell using the 

8 mirror symmetry lines Fi, F2 and F3. All sides of the sample are placed 

9 in a uniform gas environment. Two charge-carrying species are considered: 

10 oxygen vacancies, denoted by the subscript ion, and electrons, denoted by 

11 eon. 

12 The framework we propose is very broad in scope, however we specialize 

13 our study to Samarium Doped Ceria (SDC). Doped ceria is a class of mate- 

14 rials that has recently gained prominent relevance in fuel cell technology—.— 

15 We suppose that the uniform gas environment consists of a mixture of hydro- 

16 gen and water vapor and we solve the electrochemical potential and current 

17 of both charge carriers using a linear and time-independent model, which we 
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1 develop via perturbation techniques and Fourier transformation. We mainly 

2 compare our computational work to the data of Lai et al.- but we also lever- 

3 age on some results of Chueh et al.- to justify the boundary conditions. Both 

4 works study SDC-15 (15% samarium doping), hence the background dopant 

5 particles per unit volume, B, is well defined and reported in Tab. [H 

6 The surface dimensions are kept constant: the width of the metal | ceria 

7 interface (r4) is 2Wi = 3/im and the width of the gas | ceria interface (Fs) is 

8 2W2 = 5fmi. The thickness of the MIEC is set to be 2/2 = Imm, unless oth- 

9 erwise specified. Due to high electronic mobility in the metal, the thickness 

10 of the metal stripe does not affect to the calculation, and thus the thickness 

11 of the electrolyte is in effect the thickness of the cell. Hence we assume that 

12 the characteristic length scale of the sample under study is Ic = 10/im. The 

13 data mentioned above is summarized in Tab. [TJ 

14 The assumptions of the model are rather standard for MIEC. We set that 

15 the gas I metal | ceria interface, or triple-phase boundary, has a negligible 

16 contribution compared to surface reactions.— We further treat the surface 

17 chemistry as one global reaction, and do not consider diffusion of adsorbed 

18 species on the surface.— Combined with the final assumption that the metal 

19 I ceria interface is reversible to electrons, i.e., a Ohmic condition,- we are 
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only considering two steps in the electrode reaction pathway, for instance, 
surface reactions at the active site of the SDC|Gas interface and electron 
drift-diffusion from the active site to the metal current collector both along 
the SDC — gas interface and through the SDC bulk. 

We indicate the equilibrium quantities, such as electron and oxygen va- 
cancy concentration, with the superscript (0). In order to determine equilib- 
rium concentrations of charge carriers, we consider the following gas phase 
and bulk defect reactions: 

if2(gas) + 02(gas) ^ HiO^gas) (la) 
Oo ^ VS' + l02{gas) + 2e' (lb) 

where the Kroger- Vink notation is used,— i.e. V" is a vacant site in the 
crystal, e' is an electron, and Oq an oxygen site in the crystal (superscripts 
• , ' and X indicate respectively +1 charge, -1 charge and zero charge). At 
equilibrium the number of vacant sites per unit volume is c-^^^, and the num- 
ber of electrons per unit volume is CeoL At equilibrium the following two 
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quantities will be constants: 



= #^ (2a) 

PH2P02 

(0) \ 2 ^(0) 



= "-fP'£ (2b) 



in addition to that, electroneutrality will be satisfied throught the sample: 

(0) (0) 

1 + ^ - 2%^ = (3) 
B B 

where pk = and pk is the partial pressure of species k. In the di- 

latm 

lute limit, at a given temperature and partial pressure, we solve for the 
equilibrium concentrations of vacancies cf]^ ~ B/2 and electrons cioL ~ 



r> V2Kr 

Finally we assume that the mobilities u of all species are given in Tab. [21 
from Lai et al.- 
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3 Background 



3.1 Asymptotic Modeling of Mixed Conduction in the 



A mixed conductor is a substance capable of conducting two or more charged 
species of opposite sign. Mass and charge transport in sohds are described, 
at a mesoscopic level, by drift diffusion (DD) equations. The derivation of 
these equations is given in many textbooks, see for example.— For clarity 
we will shortly rewrite them here. For a mobile species m, the continuity 
portion of the DD equations is expressed by equations of the form: 



where is the concentration of species m, is the particle (superscript 
P) flux of species m per unit area and Um is its net rate of creation per unit 
volume. 

We will assume the following phenomenological relationship for the flux 
of species m (this relation is valid for VT ~ and VP ~ 0—): 



Bulk 



dc. 



m 



(4) 



dt 



m 




m. 



(5) 



hT 



m 
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where Dm and jlm are respectively its diffusivity, given by Einstein's relation 
Dm = UmkbT/zm {um is the mobility), and its electrochemical potential, 
given by an expression of the type: 

jlm = ^^m + hT\og{Cmfm{Cm, T, P)) + 2:^60 (6) 

In the latter e is the elementary charge and (p is the electric potential, fm 
is the activity of species m and Zm is its integer charge, i.e. -1 for electrons, 
+2 for oxygen vacancies in an oxide and fi^ is a reference value. We also 
define the T<r-electrochemical potential of a species m as: 



f^l = (7) 



The same equations are sometime expressed in a different way; if we define 
the conductivity am = ^-^^j^f^^ we will deduce from Eqn.sS] and [H] that: 



V ■ <[Dm+ TTj VC^ H V0 > = UJm (8) 



at l\ a log Cm J 

Here we suppose the presence of two mobile species: oxygen vacancies, 
which we indicate with the subscript ion {zion = +2), and electrons, subscript 
eon {zeon = The distribution of electrons and vacancies is thus described 
by 3 equations: one for the electric field (Poisson's equation for the potential) 
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and two for the mobile species conservation. This set of equations can be 
written as: 



A(f)= -{B + c^on - '^Cion) (9a) 

dtCeon + V ■ (^-DeonCeon'^J^^ = (9b) 
dtCion + V ■ ( -DionC^on^^^^ ) = (9c) 



where e is the permittivity of the medium, B is the background dopant 
concentration in number of particles per unit volume and where we have 
chosen Weon = 'jJion = 0. In the dilute limit 171819 ^20 j^g^g. 

jleon = ksT log ( ^ J - e0 + (10a) 

\ '-'eon / 

fi^on = kBT log (^) + 2e4> + ill^ (10b) 

1 where c°,„ and Cg^^ are reference values. 

Non-dimensionalization of the Eqn.s[9]with respect to its relevant param- 
eters proves to be crucial in order to understand appropriate time and length 
scales. We apply the transformations: (x, t) — » (x, t) such that x = /gX and 
t = rt. At this point we suppose the diffusivities -Dgon and Z)jo„ are uni- 
form (we shall relax this approximation later). Also, we define Ut = kbT/e, 
(j) = (P/Ut, Tn = ll/Deon, Tp = l^/ Dion and T = min(r„,rp). Obviously 
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V^(-) = and = ^di{-). So Eqn. [9] becomes: 

A,</) = ^ 1 + - (11a) 



ion , 



T 



Ceon Ceon 



7^*7or - ■ 1 v^'^r + 1 = (lie) 



^lon ^ton 



where cion and cf2i are equilibrium values.- Define now the Debye length 



= y and \ = j^. We suppose A 1, which holds true for highly 
doped MIECs and sufficiently large characteristic dimensions, and we use 
singular perturbation of Eqn. Illal to obtain:— 

(0) (0) 
^ Ceon Ceon _ r) ^iora ^ion _ p 

^ B 5 ~ ^ 

Leon (-jo„ 

In view of the latter, we can drop Eqn. Illa[ thus we are left with Eqn.s [TTb|lllcl 
and [121 We now focus on impedance conditions, i.e. we suppose an off- 
equilibrium perturbation of the boundary conditions which in turn will slightly 
affect all unknowns (terms with superscript (1) are much smaller than the 
terms with superscript (0)): 



(13a) 



(0) , (1) ^ (0) 1 , nouN 
eon ^ ^eon ^eon \ ~ (Q) \ \^^^) 



Ceon '^eon ' '-'eon '-'eon 



c- - c^°^ + c^'^ - c^°^ f 1 + ^ 1 (I3c) 

^lon — ^ion ^ ^ion ~ Hon \ ^ ^ (VC\ \ \^'^^) 



ion , 
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We set n^^^ = ^fgf- and p^^^ = and suppose CeJn, c-°„ are uniform and 

''eon '^ion 

(f)^^^ = 0. If we also use the definitions of Eqn. [13] in the Eqns. Illbl and lllc[ 
we obtain: 



^dt (1 + n«) + ■ (-V,(l + nW) + (1 + n«)Vi0«) = (14a) 
^dt (1 + - V, ■ (V,(l + + 2(1 + V,0(i)) = (14b) 

If we retain in Eqn. [14] only first order terms, we get: 



T 



(15a) 
(15b) 



The electroneutrality condition, Eqn. [12], at first order gives that p 



(1) 



-n 



= ——n''^\ Thus defining: 



2 c(°) 2 p 



T + — r 



1 + 



4p 



1 + ^ 



(16a) 
(16b) 
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helps rewrite the Eqn. [TS] as: 



^diu^^^ - A^n^^) = 
r 



(17a) 
(17b) 



3.2 Boundary Conditions 

It follows from symmetry, Fig.[T], that dxcf)^^'^ = dxU^^^ = on r2 and T^. Since 
the metal is ion-blocking, ^"^dyvS^^ + 2dy4'^^^ = will be satisfied on We 
assume as well that the response of the metal to an electric perturbation is 
fast compared to the MIEC, from this it follows that we can take the electric 
potential (f>^^^ uniform on Thank to linearity and given the impedance 
setting, we can choose 0*^^) = ^e*'^'^*j on r4 and (p^^^ = n^^'' = on Fi. 

We assume the chemistry due to the reactions on has a finite speed 
and that it is correctly characterized by a one-step reaction.— For simplicity 
we start from: 

i/2(gas) ^ //20(gas) + V^' + 2e' (18) 

We also remark— that the rates of injection of vacancies uJion,s and elec- 
trons u!eon,s ^t F5 Satisfy (subscript S indicates surface) the following two 
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equations: 



^ion,S o^eon,S 

(19) 

^ion,S = ^fpH2 ~ ^rpH20CionC^eon 

where kf is the forward rate of the reaction in Eqn. [TS]and is the reverse 
rate. 

The latter gives, under small perturbation assumptions,— a Chang- Jaffe 
boundary condition:— 

(1) _ . Aon r 1/4 [ ^ \ ~ (1) 



- -\l^s = ^^k^c^: [ 1 + ^ ) PH.n^'^ (20) 

where we choose 



We suppose kj = 2^kf and kf = kjp^o^ ^ ^^^S^^ 
/? = 1/4.^ 

Hence the y-fiux of electrons and vacancies satisfies the following ex- 
pression along Tg: jf„„ ■ By = 2j^„ ■ = -Ueon.s- If we define = 



^ ■f^'M ~ n^°" ) ^^"^ ^ ^f^~W ~^ 4 I ^ we can rewrite the 



boundary conditions on on Fs as dycf)^^^ = A^n*^-^^ and dyrS^^ = An 



n 



(1) 



^ The order of magnitude of kf is given by: = 10 "^"^^ h ~ 10 ^m, I?iori = 

10-lOmVs and fc? « 10^2, so fc. « 1032 X ^ X 10-6 = 1021 #par»clc. ^ 10-3 mo;. 



j^n — 7 mol 
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1 3.3 Weak Formulation of the Model 



If we Fourier transform Eqn.s [T7] and the boundary conditions with respect 



to t 



we find the 



formed quantity) If 



'ollowing system of equations ((■) indicates Fourier trans- 
which we call IS equations: 



*n^^) _ A<^(^) = 



(21a) 
(21b) 



on Fi 

on F2 A F3 



with boundary conditions: 

0(1) =0 A ri(i) = 

9^0(1) =0 A S^n^i) = 

0(1) = 1 

a^0(i) = i0n(i) A dyfi^^^ = i„n(i) on F5 
We can recast the Eqn. [21] and [22] in weak form taking as test functions 



A 9^n(i) = -4f 9^0(1) on F4 



(22) 



^We choose unitary Fourier transform f{uj)~ f{x)e '"^ dx 

^We factored out the Dirac distribution that comes out of Fourier transformation of an 



exponential 
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3 m^e, eH\n\ Ti), z^,^ G ff^^ \ (ri U r4)):^ 



^< I n^^i^Re dA - Vn^^l ■ VniRe dA+ i^n^J^mRe dA + 



rs 



n 



(23a) 



r4 



/ nj^e^imdA + / Vn^^j, ■ Vm/m - / Anfijlmimdx 
Jn Jn Jrs (23b) 

+ 4^ / aj^05'im,„dx = O 

/ ^imV'iJe di - / V0g^ ■ VijRe dA + [ i^n^j^ije dx = p3c) 
Jn Jn 

I n^T/;,^ di + / V0i'i ■ V^,™ di - / i^ni^i^A.^ dx = 023d) 
Jn Jn Jv^ 



with the condition that: 

0^=0 A (^«=0 on V, (24a) 

ng^ =0 A 411 = on Ti (24b) 

A 0gi = O on 14 (24c) 



= 1 ^ ^^^^ 



2 It is easy to show that the sum of the Eqn.s [2S1 is bounded and thus the 

3 bihnear form associated to the weak formulation of Eqn.s [21] with [22] is con- 

4 tinuous. Further, the problem is weakly-coercive hence it admits one unique 

5 solution.— 
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3.4 Numerical Solution Procedure for the 2D Case 

In order to solve numerically the Eqn.s[23]with boundary conditions Eqn.slMl 
we employ an h-adapted finite element method (FEM), implemented with 
FreeFem++.-^ The governing equations are discretized on a triangular un- 
structured mesh using quadratic continuous basis functions with a centered 
third order bubble. We use a direct method to solve the linear system fol- 
lowing integration of Eqn.s [23] in the discretized mesh. Then the mesh is 
adaptively refined nine times for each case. The a posteriori adaptation is per- 
formed the first six times against the 4 dimensional vector ^V3fJ 
and subsequently against r^^, see Appendix |X1 The h- adaptation ensures high 
regularity of the a posteriori estimator,— locally below 10~^, and it guar- 
antees that the mesh is finer where sharper gradients occur. Independently 
of frequency, mesh adaptivity results in coarseness everywhere except in the 
vicinity of the interfaces, in particular the refinement increases towards the 
triple-phase boundary (the intersection of metal, oxide and gas phases, which 
is though to be a particularly active site for electrochemical reaction s -"^^ ^^ ); 
this fact indicates strong non-linearities around that area. Finally we note 
that FreeFem-|-+ execution time is comparable to custom-written C++ code 
and its speed is enhanced by the utilization of fast sparse linear solvers such as 
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fJ^eon , V JL ftjo„ 



12 the multi-frontal package UMFPACK.— Due to the sparsity of the problem 

13 we make extensive use of this last feature. 

14 We further note that the utilization of asymptotic expansion and Fourier 

15 transformation techniques, while guaranteeing linearity, has a great speed 

16 advantage over direct sinusoidal— and step relaxation techniques.— Further, 

17 this method can be directly used to examine chemical reactions within the 

18 cell and draw directly conclusions about fast and rate-limiting chemical re- 

19 actions. Also, this procedure lends itself to direct error estimation and its 

1 implementation can be done automatically for a time-dependent problem.— 

2 3.5 ID case: Analytical Solution 

3 Since we also aim at comparing the ID and 2D solutions, it is beneficial to 

4 revisit the ID solution of Eqn.s[2Il- The solution {h^^\ c/)^^'') will satisfy (if 

5 Uy^O): 



n 



a±e 



(25a) 
(25b) 
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where for simplicity we indicate = e*4. The boundary conditions, as in 
the 2D case, at y = (Fi) are: 

0(1) = A n(^) = (26) 

The latter can help rewrite Eqn.s [25] as: 

n*-^-* = 2a_|_ sinh (^y/iy^r^uy^ (27a) 



y + 2a+^ sinh fViy^y) (27b) 



If we set 70 = ^^^-^'^i^y and 7„ = i|70, then at y = k we have the 
following conditions:- 



^W = l A nW+70^ + 7n^ = O (28) 

The boundary conditions Eqn.[2S]will lead to the determination of a+ and 
^0^^)' in Eqn. [27] and the ID model leads to impedance of the form-^:— 



, - ■'jon ~l~ Reon . , 

s H —J tanh s 



(29) 

where all the relevant terms are reported in Table [3] 



21 



4 Results 



6 4.1 Comparison to Experiments 

7 The electron electrochemical potential drop across the sample, i.e. the elec- 

8 tron electrochemical potential difference between the top and bottom elec- 
1 trodes and its symmetric reffection), is given by the following expression: 



(30) 



where < a >\ indicates the average of the quantity a over the set A. At first 
order the ^-electrochemical potential is given by (^f^'e^^ = 4>^^^ ~ fi'^^^- The 
electric current density at the the two ends of the circuit is: 



ill) 



{W^ + W2)lc 



Hence, the 2D impedance is given by the expression: 



(31) 



(32) 



We define the error of the 2D impedance Z2D with respect to experimental 
impedance Zid spectra Eqn. [22] as follows: 



1 



Z2d{uj,po2,T) 



Zir,iLU,po2,T) 



(33) 
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For every data point, uniquely defined by the couple {po2,T), we fit the 
2D data against the measured ID equivalent circuit data in- by minimizing 
eF{uj,po2,T) with respect to the surface reaction constant = ApQ^, which 
is a function of both02 partial pressure and temperature. We remark that 
is the sole parameter we allow to vary in this procedure and all other 
data is obtained from the literature and presented in Tab. HI With only one 
parameter variation, we obtained excellent agreement between experimental 
results and 2D calculations, i.e. eF{i^,Po2^T) < 2%. As an example, 2D 
results at four different oxygen partial pressures and at 650°C are shown 
in Fig. [3l We computed the /cj by minimizing the ep for a total of 28 
cases (7 pressures times 4 temperature). We report in Tab. |4]the results of 
linear regression of these minimizing values (each line is derived on keeping 
the temperature fixed and varying ^02)- We also write in Tab. HI the 95% 
confidence intervals for the fitting of A, i.e., A ^ A ± ea, and a, i.e., a = 
a ± Sa', we finally report the root mean square error a and the adjusted R- 
squared,— regarding the latter, a value close to unity indicates a perfect fit 
while negative values indicate poor data correlation. Directly from analysis 
of Tab. m we deduce that k^ fitting to a straight line is reasonable for "high" 
temperatures (T > 550°C). We note that k^ is temperature-dependent via 
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A {A decreases with T). Furthermore k'j is shghtly pressure dependent via 
the coefficient a, the average value of a ~ 0.05 > 0, however the error is of 
the same order of the slope. Hence the total rate of reaction is very likely 
to be UiJeon,s Pol^^~^^ where /? is somewhere in the set [0,0.1], most likely 
equal to 0.05. 

4.2 The Polarization Resistance in Frequency Space 

One of the goals of fuel cell science is to understand and possibly reduce the 
polarization resistance, i.e. that portion of the resistance due to electric field 
effects at interfaces. For that purpose it is key to identify and understand the 
main processes that intervene in the definition of this quantity. Specifically, 
the area specific polarization resistance for our system is defined as:— 

y± TT < Aion >r5 - < Aeon. /o/|\ 
^ion - TT^N (,<54j 

JlP 

where j}p = /rg ^eon,s is the ionic contribution to the area specific 

current. The can be understood as the sum of a surface Zg^rf and a bulk 
polarization resistance, Z^uik = Z^^ — Zgurf, where the Zgurf is the portion 
of the area-specific resistance due to effects of the exposed boundary and 
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it is given by: 



'surf — Ut 



< f4on >Ts - < A, 



'eon 



(35) 




In our model, by definition, the Zg^rf ^ IR"^ is proportional to (1 + W1/W2) 
and inversely proportional to both and kj: 



impedance is due to surface effects. From Fig. H] we note two fundamental 
facts: first, as we expect, at "lower" injection rates the fsurf increases, physi- 
cally this means that that if the chemistry is sufficiently slow it will dominate 
the polarization resistance leading to an fsurf of approximately unity. Sec- 
ond, we notice frequency dependent behavior of -R^„. Our computations 
show that /surf decreases with u, while the dephasing between Z^urf and Z^^^, 
described by arg(/surf), increases with /cj and decreases with uj. The behav- 
ior of /surf in phase space clearly shows that Zsurf includes two interrelated 
processes: 

1. reactions on the surface exposed to the gas; 

2. transport of charged species in MIEC. 




(36) 



The fraction fsurf 



surf 



indicates what portion of the polarization 
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14 Within this framework, as uj increases, the losses in the polarization due to 

1 drift diffusion increase and surpass the (constant) reaction or surface losses. 

2 4.3 Analysis of the 2D Solution 

3 4.3.1 Qualitative Considerations 

4 We can then use the framework to study the two complex electrochemical po- 

5 tentials /teon = n^^^ — ^^^^ and /tion = (f)^^^ + ^n^^^ as functions of frequency. In 

6 Figs. [5] and [6] we plot the 2D distributions of the latter in the computational 

7 domain at T = 650°C, = 10^^^ and = 10^^ with frequency u increas- 

8 ing from 10"'^ to 10^ rad/s . Thank to the Figs. [5] and El we can address the 

9 qualitative behavior of the solution. We first analyze the qualitative distri- 

10 bution of fluxes: from the gradient of l/ieonl, which gives an idea of electron 

11 flux, that electrons flow from the gas|ceria interface Fs onto the ceria|metal 

12 interface F4 through a cross-plane current Ig'^ , and concurrently electrons 

13 flow onto the ceria|metal interface F5 from its mirror symmetric counterpart. 

14 Similarly the M1EC| metal interface is blocking to vacancies, hereby the va- 

15 cancies correctly flow from the bottom to the top ceria|gas interface F5. It 

16 is also clear that the complex potential of the electrons /ieon changes signif- 

17 icantly as u increases, while /iion is relatively unaffected. The penetration 
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depth, which is defined as the vertical displacement from r4 where surface 
electrons can penetrate into the bulk, decreases with uj as the ID model 
hints (in Eqn.s [271 the solution decays exponentially with 1/ ^t*uj). As uj 
increases, the dephasing of /teon first increases and then decreases and it is 
weakly dependent upon the distance from or conversely, the penetration 
depth into the MIEC. We notice that the same dephasing increases and then 
decreases for fiion- However, while for the vacancies, the behavior of \(iion\ 
and aig{jlion) is qualitatively the same, this is not the case for the electrons, 
where through a wide array of cu's, the qualitative behavior of \jleon\ and 
arg(/ieo„) is distinctly different. 

Deriving the electronic and ionic currents from the computations requires 
some care and it will not simply be V|/im|For example, for electrons, we note 
that: 




(37) 



We will call the complex current j' 



;C 



J 



= c, 



(38) 



eon 
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the physical current will be l 



J eon 



^ OL) 



(39) 



In order to compare the ID and 2D solutions qualitatively, we first focus 
on the case on = where = 10^^, and we shrink the size of the slab while 
keeping the same framework and model parameters. This corresponds to a 



in Fig. [H] the results of the computations in the case where the conditions 

are very reducing. We depict what happens to Rion, Reon, Rton ^"^^ fsurf as 

AR changes. We notice that decreasing AR corresponds to an increase in 

effective electronic and ionic resistance compared to the ideal case computed 

according to Tab. [3] which in turn corresponds to AR oo. Deviations 

from ideality occur already for AR ^ 25, hence even for reasonably large 

AR the ionic and electronic resistances deviate from the ideal ID case, this 

is clearly shown in Fig.s [8] a and b. The same applies to the polarization 

resistance -R^„, Fig. [Hb, which is fiat above AR ^ 25, below this value 

sharply increases due to bulk polarization effects. As the deviation from the 

ID setting starts, not only ionic and electronic resistivities change, but so 
^We remark that for complex valued function fj. in general we have abs(V/i) ^ 



decrease of the aspect ratio of the sample defined as AR 



. We show 



W1+W2 



V (abs (^)) 
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does the relative importance of surface and drift diffusion effects. Hence the 
polarization resistance is thickness-dependent, and the dependence is due 
to the emergence of two-dimensional effects. The increase in drift diffusion 
resistance due to the motion of electrons from Fs to r4 is also shown in 
the fsurf which increases with the AR reaching unity for AR — * oo. This 
effect is even clearer if we plot the electrochemical potentials of electrons and 
vacancies at u; = 0, we note a shrinking of the affected the sample 

thickness decreases corresponding to an increase of polarization resistance. 
This effect is purely 2D and cannot be studied using a ID model. 

4.3.2 Quantitative Analysis 

In order to compare the ID and 2D solution quantitatively we define the 
following two functionals: 

^ j \fiiD{y',uj) - fi2D{x,y',Lj)\ dx 



C[fllD,P'2D,y,Uj] 



{fiiD{y\ u}) - fi2D{x, y', uj)) dx 

y'=y 



{40b) 



The functional v describes the "pointwise" distance between ID and 2D 
solutions of /i at a section y and the functional C describes the "average" 
distance between ID and 2D descriptions. Physically v indicates how far 
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11 apart the ID and 2D electrochemical potential are, while ( "measures" the 

12 soundness of fitting a ID case with the 2D model. We can examine the 

13 applicability of the ID approximation for data fitting via (. 

1 In order to further compare the 2D model and ID model and demonstrate 

2 the importance of 2D effects adjacent to the injection sites, the "pointwise" 

3 distance u and the "average" distance ( defined by Eqn.s liObl are computed 

4 at the same conditions (T, poj, k^) in the frequency range of 10^'^ ^ uj < 10^ 

5 rad/s along the symmetry axis T2, Fig. [91 In the first line we plot the case 

6 where the sample is very thick with respect to the horizontal dimension 

7 {AR = 125), both the Ueoniy,^) = ^ [f^eon,iD, f^eon,2D,y] and the (ion{y,uj) = 
3 C [f^ion,iD, fJ'ion,2D,y] are cxtrcmcly small and the adjacency between ID and 

9 2D impedance is near perfect. If we decrease AR to 12.5, then the ID and 

10 2D solutions tend to be further apart with Ug ^ 25% and Ce up to 20%. The 

11 difference between the two further increases at AR = 5 where the difference 

12 between impedance spectra is significant. 
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5 The Effect of Diffusivity Gradients 



5.1 Extension of the Model 

Interface effects are one of the biggest sources of uncertainty in doped ion- 
ics because impurities in doped materials tend to segregate near interfaces 
and affect electro-catalytic processes, absorption and diffusivities near the af- 
fected interfaces. Many studies have attempted to address these issues. 
However, to the authors' knowledge, no continuum model has addressed yet 
the relationship of these changes to polarization resistance nor to impedance 
spectra. In this part of the paper we intend to address the effects of non 
uniform diffusivities, which are localized near the interfaces, and which we 
imagine are due to impurity segregation at the exposed surface in Fig. [1]) 
and to the MIEC|metal interface (r4) 

We shall assume that diffusivities near the MIEC|Gas interface and MIEC| Metal 
interfaces have non- zero derivatives only along the y direction. We further as- 
sume that diffusive effects are symmetric on both ends of the sample y = ±^2, 
hence do not affect our initial symmetry assumptions. Lastly we suppose that 
the functional form of the diffusivities are known in the MIEC and are given 
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by: 



Dl = ^+ yhlK - IJ (41) 



where m can be either eon or ion, and A^,, the length scale of diffusive 
changes, is much smaller than Zc, the characteristic length-scale of the sample 
(Am ^ /c)- We stress again that the main assumptions are that the diffusivity 
gradients parallel to the interfaces are null and that the diffusivity gradients 
do not affect bulk properties of the material nor the defect chemistry. In 
other words, near-interface effects involve only diffusivities. 

Under the same small perturbation assumptions we used above we can 
deduce that the equations that describe the impedance spectra behavior of 
ions and electrons are given by □: 

n(^) = ^p(i) (42a) 
P 

+ V, ■ (-Z^:„„ (V,n(^) - V,0(^))) = (42b) 

d-tP^'^ + V, ■ (V,p(^) + 2Vi0(^))) = (42c) 



The sum of the Eqns. |l2] and their weighted difference lead to (see appendix 



In order to ensure linearity, we assume that Dfen(i'V0(i) < IZJfcVn^i)! w ZJfcV^ 
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ID: 



where: 



^din^^^ + ■ (-aiiViU^^^ - aiaV^^^^^) = (43a) 
diu^^^ + Vj. ■ (-a2iVin(i) - 022 V^^^^^) = (43b) 



T 



1^ ~ 1 I A ' ^2 TTA ^ ^ 

- D* D* +^D* 

«21 - — 1| , ^22 - — ^ (44bj 

n fi 

2 The Eqn.s |13] with appropriate boundary conditions, Eqn.s [221 are quasi- 
1 hnear and hence can be Fourier transformed. In short they can be recast in 
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2 weak form as in Eqns. [221 



/ '^fi^Re dA - an VnJJ] ■ VmRe dA - ai2V0^^] ■ Vm/je dA + . . . 
Jn Jn Jn (45a) 

• • • + J An,2n^^rnRe - 4^ j dyctl-^^mRe dx = 

^'^n I fi'Rl'f^ini dA + O'u'^nfi " Vm/^ dA + ai2 V0^^j, ■ Vm/^ dA - . . . 

7n 7n (45b) 

• • • - / An,2nfimi^ dx + 4- / dy(f)flmj^ dx = 
^r;j^4lijRedA - ^a2lVf^g^ V^Redi-^a22V0g^VV^Redi 

(45c) 

+ /r^ A^,2n'ReijRe dx = 

^^;L^S^/-di + /f,a2iV4;i-V^,^di + /^a22V4;i-V^,^di 

(45d) 

- /r^ ^</.,2^?i^7m dx = 

3 

where: 

^n,2 = a-iiAn + ai2^0 (46) 

^^,2 = a2lAn + 012^0 (47) 

(48) 

4 If we change the diffusivity of vacancies at the gas|ceria (Fs) and metal|ceria 

5 (r4) interface by changing aion, we need to adjust the as follows, in order 
1 to keep the same rate of injection Cjf^, Eqn. [2D1 
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~kf = ™^ (~kj) (49) 

2 Numerically we use the same approach described for the linear case but we 

3 need the error estimator to account for off-diagonal and space dependent 

4 parameters, Eqn.slH] (in the linear case an = 022 = 1, ^12 = ^21 = 0). 

5 Finally we note that we assume that the model holds for length-scales 

6 just one order of magnitude greater that the lattice parameter.— This ap- 

7 proximation can be justified heuristically using the work of Armstrong—,— 

8 which shows that deviations of the continuum drift- diffusion approach from 

9 atomistic models are usually small, even in cases where field effects are big. 



10 5.2 Results of the Model 

11 We first ran the model at steady state {u = 0) with the objective to an- 

12 alyze the fsurf = at u; = for a wide array of parameters aeon = 

ion 

P)SURF / r^BULK j _ r^SURF / p)BULK „,V,„^„ ^ _ ^ „„rl \ — 

13 Ueon I ^ eon aiiU. tXion ^ion l^ion i Wlieie Clgon u^ion aiiU. neon 

14 ^ion ^t varying kj^\ For reasonable fitted values (Tab. Hj) and for a wide 

15 parameter set, we show that the polarization resistance is surface dominated 

16 making fsurf ~ 1 robustly. 

1 If chemical reaction rates are "sufficiently" slow (e.g. /cj ~ 10'^^) and if 
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2 the sample is sufficiently thick, then the polarization resistance is dominated 

3 by surface effects in the linear case (ajon = 1)? corresponding to an absence 

4 of diffusive gradients at the exposed surface. If impurities are present at the 

5 exposed surface, diffusivities of charged species may change and hence one 

6 could argue that the polarization resistance is not surface-dominated. In or- 

7 der to address this point, we ran two limiting cases, one featuring "slow" 

8 chemistry (/cj(ajon = 1) ~ 10^^) and the other one at "fast" chemistry 

9 {k'f^aion = 1) ~ 10'^^). We present the results of these calculations in Fig. [10] 

10 where we plot fsurf as a function of both aion = o^eon and the diffusive gra- 

11 dients Xion = Kon- We notice from Fig. [TOk that fsurf is very close to unity 

12 for two order of variation of surface-to-bulk diffusivity ratio 0.1 < a^on ^ 10 

13 and for a wide span of diffusivity length-scales 5nm < A^on < 1/im. This 

14 indicates that if we perturb the the surface diffusivity up to one orders of 

15 magnitude with respect to its bulk value its impact on polarization resistance 

16 is minimal. The qualitative effect on the impedance is also small as shown 

17 for a variety of cases in Fig. [11] 

18 If we choose a "fast" chemistry condition instead, e.g. /cj ~ 10"^^, the 

19 situation changes significantly from the base case {aion = 1), Fig- [TTb . In 
1 this Figure we focus on points A through D. (Pt. A), having aion = 0.1 
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2 and Xion = 5nm, indicates that near surface diffusivities are an order of 

3 magnitude lower than their bulk value and this deviation is concentrated 

4 near the surface: in this case the polarization resistance is drift-diffusion 

5 dominated. If the diffusive length scale is increased to Xion = ^f^'Tn, while 

6 keeping ctjon = 0.1, (Pt. B), the fsurf will not decrease much further. Starting 

7 from (Pt. A) we can move to (Pt. C), where diffusivity gradients are sharp 

8 {Xion = 5nm) but the diffusivities at the surface are an order of magnitude 

9 greater than its bulk value. In this case, the fsurf increases because of the 

10 increase in the bulk diffusivity. Going to (Pt. C) to (Pt. D) increases the 

11 length-scale of the diffusive effects leading in turn to bigger increase of fsurf- 



12 We can summarize our findings as follows: 

13 1. if the rate of injection of electrons is sufficiently "small" (slow chem- 

14 istry) and of the order of the fitted values reported in Tab. HI then 

15 the diffusivity grandients localized at interfaces will affect little the 

16 polarization resistance and the impedance spectra; 

17 2. if the chemistry is sufficiently fast, sharp changes in diffusivity can affect 

18 strongly not only the impedance behavior but also the polatization, 

19 in particular if the diffusivities increase sufficiently, strictly near the 
1 interfaces, the polarization effects will shift to be surface dominated, 
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2 while a decrease is associated to drift- diffusion dominated polarization 

3 resistance. 

4 6 Concluding Remarks 

5 A general two-dimensional numerical framework has been developed for the 

6 coupled surface chemistry, electrochemistry and transport processes in mixed 

7 conductors based on the finite element method. As a specialized application 

8 of the framework, a time- dependent model was formulated based on first- 

9 principles for the AC impedance spectra (IS) of a samaria doped ceria (SDC) 

10 electrolyte with symmetric metal patterns on both sides, and the IS was sim- 

11 ulated for typical fuel cell operation conditions in a uniform gas atmosphere 

12 {H21 H2O, Ar) at thermodynamic equilibrium using the small perturbation 

13 technique. 

14 The validity of the model is demonstrated by fitting to experimental (ID) 

15 impedance spectra data of an SDC cell in literature, varying only the reac- 

16 tion rate at the SDC-gas interface. Excellent agreement (< 2% error) was 

17 obtained. We then numericalUy investigated the influence of the variation of 
1 several parameters on the polarization resistance and the impedance spectra, 
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2 especially within regimes not probable for the ID studies. Our calculation 

3 shows that the 2D effect of cell thickness variation on the spectra becomes 

4 pronounced as the aspect ratio goes below a certain threshold (25 for this 

5 work); surface reaction dominates the polarization resistance when the in- 

6 jection rate at the SDC surface exposed to gas is sufficiently slow; sharp 

7 gradients in diffusion coefficient strongly influence both impedance behavior 

8 and polarization when surface chemistry is sufficiently fast. 

9 The discussions in this work provide useful insights into the correlation be- 

10 tween materials properties of SDC and its applications in fuel cells, intensely 

11 studied by the solid oxide fuel cell researchers. In addition, the geometric 

12 capability (up to 3D) and high computation efficiency makes this numerical 

13 framework an ideal tool for the general study of mixed conductors. 
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2 A Error Estimator and Refinement Strategy 

3 The local residual for nRei at a triangular element K of the mesh, can be 

4 computed as follows:— 



K 



e,hl ''^'^n'^Im,h 



ail 



+ ai2- 



(1) 





j 




JVir\K 



I 



dn 

(50) 



dn 



(r2ur3)nx 



"'K 



5 where |a] is the jump of the quantity a across the faces of h is a. 

6 measure of the size while hx is the measure of the size of the sides of 

7 K. Similar residuals can be found for n^^, 0^^, 0^^. Their sum ^^.r^ 

8 constitutes a reasonable local a posteriori error estimator. is a weakly 

9 coercive upper bound for a||M||L2 ~" ^11 VmH^j where a and h are constants and 

10 a — \ Ai^g, Aij-^, (/^j^ 



11 B Derivation of the Non-linear impedance Spec- 

12 tra Equations 

13 We start with the electro-neutral form of the drift-diffusion equations, where 
1 we assume that the diffusion coefficients normalized with respect to their 
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bulk value = D^^urf/^bulk. 

^a^-^d) + V, . ^-D:^ (V,n(i) - V,0(i)) ) = (51a) 
+ V. ■ (-Z^L V.n(^) + V.0(^)) ) = (51b) 

We first sum the Eqn.s ISTal and [STbl and obtain: 

(52) 

Finally we multiply Eqn. I51bl by ^ and sum to Eqn I51a[ 

- 7) 5*-«+V.- (- (DL - DU V.r^(^) - (d:„„ + V.0«) 

(53) 

From the and , the Eqn .sHSl follow immediately and so do their coefficients 
given in Eqn.s HH . 
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Table 1: Data for the domain geometry and background doping 
Wi 1.5 fim 

W2 2.5 fim 

I2 500 iim 

Ic 10 /im 

B 3.47 X 10+27 #particles 



Table 2: Temperature range and material constants for the simulations. 



T 




550°C 


600°C 


650"^ 


K, 


5.059E+27 


4.814E+25 


7.757E+23 


1.944E+22 


Kr 


5.008E-22 


2.263E-20 


6.610E-19 


1.340E-17 


^eon 


V^s 




4.762E - 8 


6.257E-8 


6.873E-8 


8.123E-8 




2 

m 




1.166E-9 


2.070E-9 


3.359E-9 


4.936E-9 
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Table 3: Definitions of tlie terms in tlie ID model 



ion 


Measured 


ID 


Z I ^ eon 


Rion 


2^2/ Cjon 


Ro 


1/ (l/i?eon + 1/ {Rion^IRi^) 


Roo 


1/ (l/-Reon + l/-Rjon) 


Cchem 


^2/2/ (^l/(^eonCiol) + l/(^j^onc!oi)) 


D 


(.^Rion ~l~ Reon) Cchem) 


s 





Table 4: Fitted values of k^r = ApQ , 95% confidence interval 





logio ^ 


logio 


a 


£ct 




cr 


500 


32.48 


0.150 


0.05349 


0.1655 


-0.0439 


0.1577 


550 


32.10 


0.045 


0.04160 


0.0482 


0.7622 


0.04589 


600 


32.02 


0.055 


0.06674 


0.0637 


0.5378 


0.06067 


650 


31.95 


0.055 


0.05596 


0.0623 


0.4981 


0.05938 
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Figure 1: Schematic depiction of the domain under study with annotation of 
the boundary names and dimensions. The domain is composed by an MIEC 
slab of half-thickness I2 which is mirror symmetric with respect to Fi. On top 
of the slab there is a metal stripe infinitely long deposited over the surface 
F4, the surface F5 is exposed to the gas phase. The overall sample is mirror 
symmetric with respect to F2 and F3. 
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Figure 2: Depiction of the currents in the MIEC. The superscript CP in- 
dicates cross-plane current and the superscript IP means in-plane currents. 
The subscript g indicates that the the flux is due to electrochemical reactions 
at the gas|ceria interface, while the subscript e is for electrode to electrode 
current. We notice we will have four currents: one, the cross-plane electron 
flux I^'^ from the bottom to the top electrode, two the cross-plane ionic flux 
from top to bottom gas|ceria interface Ig^ and the in-plane electronic fluxes 
Ig^ from the gas|ceria interfaces to the electrodes. 
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Figure 3: The triangle indicated fitted computations while the solid line is 
the experimental value. The results are presented at 650°C varying the po^ 
partial pressure. 
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Figure 4: Plot of fsurf = ^T^ as a complex function of u. We present 

'ion 

two cases, both at 650°C, the one to the left at very reducing conditions 
P02 = 10~^^'^^ and the one to the right at P02 = 10~^°'^^, parametrized 

r(o) 

versus k) . 
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Figure 5: Plots of the complex electrochemical potential of electrons 
fi'eon{x,y,uj) as a function of x and y in the case where T = 650°C and 
P02 = 10"^^. In the top panels we depict its absolute value \fieon\ while 
at the bottom we show its argument arg(/ieon.)- The applied frequency is 
increased from left to right, going from 0.001 rad/s to 1 rad/s. 
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Figure 6: Similar to Fig.[5l we depict the complex electrochemical potential of 
ions fiionix, y, u) where at the top we show \flion\ and at the bottom arg(/ijo„). 
The conditions are the same as Fig. O and so is the frequency range. 
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Figure 7: Potentials and current lines under small bias excitation, i.e. 
impedance at = 0, at T = 650°C and = lO"^^''^'^. The fieon (left col- 
umn) and fly (right column) along with their current lines are plotted. Each 
row corresponds to a different thickness. As I2 decreases (from top to bot- 
tom row) the area affected by surface reactions thins out; this phenomenon 
relates to an increase of the polarization resistance. 
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Figure 8: Deviation of the 2D model from ID behavior as a function of the 

aspect ratio AR = {Wi + W2) /h- We consider the case where k^p = 10^^, 

T = eSO^C and we set po^ = 10^^^-^^ip=\ow), po^ = IQ-^^-s^ (p=med), 

P02 = 10^^°'^^. The and the monotonically approach their ID 

(ideal) value if AR is sufficiently large. -R^„ increases with decreasing the 

AR while the fsurf decreases, which indicates that if the thickness is reduced 

enough, the R^^ is not just surface dominated. 
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Figure 9: Plots of the v and C's of the electrochemical potential of electrons 
(plots are shown up to 5/im from r4 and Fs) as function of y and lo and 
of the impedance spectra as the aspect ratio changes (each line corresponds 
to a different aspect ratio, 2/2 = lOOO/im, 2/2 = 100/im and 2/2 = 40/im 
correspond respectively to AR = 125, AR = 12.5 and AR = 5). A decrease 
of the aspect ratio corresponds to an increase of both u and ( and an increase 
between the (ideal) ID impedance and the 2D impedance spectra. 
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Figure 10: Depiction of fsurf in the case T = 650°C and P02 = 10" 
as a function of the ratio between near interface and bulk diffusivity, ctjon = 
Df^^'^/Dgli^^ and a,^ = Df^^^^ / D^^J^^ {a^^n = «eon), and length scale of 
the diffusive gradient A^on = Aeo„, for k'^p = 10^^ (left panel) and = 10^^ 
(right panel). 
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Figure 11: Impedance of the sample under the conditions: k^^^ = 10^^, = 
10-25.33 g^j^^ 2^ _ Q50"C, where Oeon = Oion {am = Df^^^/D^^^) and 
-^ion = Kon- The sohd hne represents the case where ajon = 1, the triangles 
and the squares indicare respectively Ajo„ = 5nm and Ajo„ = l/xm. Each 
panel corresponds to a different value of ttjon- Only small deviations occur 
from the case aion = 1 
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